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In quantum many-body theory, no generic microscopic principle at the origin of complex dynamics 
is known. In this Letter, we present a method for solving the dynamics of a rather general class 
of quantum integrable and nearly integrable systems, which leads to a natural distinction between 
regular and irregular behaviors. The basic idea is, that the dynamics of the integrable quantum 
system is described by some underlying classical one, which can be analyzed using the powerful 
tools of classical theory of motion. We apply the approach to the Dicke Hamiltonian subjected to a 
periodically driven detuning. This represents generic example of a dynamically perturbed integrable 
model. We show that scattering in the classical phase space can drive the quantum model close 
to thermal equilibrium. Interestingly, this happens in the fully quantum regime, where physical 
observables do not show any dynamic chaotic behavior. 



Introduction - Many aspects of the transition from reg- 
ular dynamics of an integrable system to erratic behavior 
of a complex system are understood in classical mechan- 
ics. On the one hand, there is the Kolmogorov-Arnold- 
Moser (KAM) theorem [IH3], that proves the stability of 
weakly perturbed integrable systems. On the other hand, 
a variety of mechanisms leading to chaos and eventu- 
ally to the ergodic exploration of phase space have been 
found. For quantum systems, there exist semiclassical 
[4] and phenomenological (e.g. random matrices [4j [5]) 
descriptions of quantum chaotic phenomena, but there 
remains an important conceptual gap between regular 
and complex behaviors. In this Letter we approach non- 
trivial quantum systems away from the semiclassical limit 
and gain microscopic insight into the emergence of irreg- 
ularity by developing a theory of dynamical perturbation 
from integrability. 

A good starting point to approach regular dynamics of 
non-trivial quantum systems are Bethe ansatz (BA) inte- 
grable models, which possess a complete set of integrals 
of motion. The exact solutions of time-independent BA 
many-body solvable systems played a crucial role in the 
understanding of various fundamental phenomena and 
concepts in physics. Famous examples are the solutions 
for the Ising model, the Heisenberg spin chain, the one- 
dimensional Hubbard model or the Lieb-Liniger gas jg]. 
However, the non-equilibrium dynamics of these models 
are rich [7J [8] and much more difficult to be calculated 
within the BA than the static properties. Formulating 
a dynamical deviation from integrability is therefore not 
only a conceptual, but also a technical challenge. 

The main finding presented in this Letter is that there 
exists an exact description of the quantum integrable 
model in terms of a classical many-body interacting sys- 
tem. Dynamical deviation from integrability for the 
quantum system can then be understood in terms of the 
classical system, for which powerful tools such as the 
KAM theorem are available. We will first derive these 
underlying classical systems for Gaudin models and then, 
as an example, analyze their behavior in a periodically 
driven Dicke model. 



Dynamical ansatz and separation of variables - Within 
the algebraic BA |9], the transfer matrix generates inte- 
grals of motion, which will be sensitive to the breaking 
of integrability. All of them depend on rapidity parame- 
ters, A rn (m = 1, . . . , M, where M is the number of exci- 
tations or particles in the system), which, in the time- 
independent case, are determined by the Bethe equa- 
tions. Our approach is to introduce time-dependent ra- 
pidities, which describe the dynamics of a certain class 
of BA integrable models with arbitrary time-dependence 
of parameters. These rapidities are equivalent to posi- 
tion variables of a classical many-body system. The first 
step in the derivation of the system of classical equations 
of motions is a separation of variables for the quantum 
mechanical wave function [3 . 

To be specific, we restrict the following discussion to 
the broad class of so-called Gaudin models [6j[T2], but the 
separation of variables can be applied to any integrable 
model, such as e.g. quantum spin chains [5 j. Gaudin 
models are relevant to a number of physical systems. For 
example, Dicke models describe cavity QED, Richard- 
son models have been applied to mesoscopic supercon- 
ductivity and central spin models to quantum dots and 
NV centers in diamonds. Gaudin- type models represent 
a generic situation in the following sense: BA-solvable 
models may depend on a number of parameters, one 
of them controls the so-called " quasiclassical" behavior. 
One can expand a generic integrable model in a Taylor 
series of this parameter. The first nontrivial term of this 
expansion is a generic Gaudin-type model. Therefore, the 
quantum phenomenon we observe here for a small but fi- 
nite value of this parameter can be expected to persist in 
general case. 

For a Gaudin model with N sites, resp. inhomogeneity 
parameters, Sklyanin [3] introduced separated variables 
Vj (j = 1, . . . , N) as operator zeros of an off-diagonal ele- 
ment of the L-matrix (see Supplementary Material). For 
a wavefunction in a representation of one can then 
introduce momentum variables Vj = (d/duj) — A(uj), 
where the concrete form of the "gauge potential" A(uj) 
depends on the symmetry and the concrete realization of 
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the Gaudin model. Generically, it has a singular struc- body potentials: 

ture, A(u) ~ Y^j=oh( u ~ z j)~ X w ^ n some given local 
spin quantum numbers lj , provided u is sufficiently close 
to the vortices (or monopoles) situated at points Zj in the 
complex plain. This structure is related to effective 2D 
Coulomb gas, similar to the situation in the Hall effect. 
In addition to these separated variables, there exists as 
well a special overall symmetry operator uq. 

The main advantage of using separated variables is 
that the transfer T-operator of the Gaudin algebra, which 
represents a generating function for the Hamiltonian and 
the conserved quantities, is written as a combination 
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of separated kinetic terms, T{uj) 

the Hamiltonian has the form of particles moving in a 
monopole gauge potential given by A{u) and the time- 
dependent Schrodinger equation for the wave function 
^({uj}) can be separated, 



iip(u) 



d 2 iP(u) 
du 2 



2A(^)^W + (A 2 (^) - A\u))^(u) 



du 



(i) 



for every Uj = with ^({uj}) = E[jLo ^( u j)- 

Classical system - Our next conceptual step is to show 
that each differential equation ( lj for one individual sep- 
uj is equivalent to a Hamiltonian sys- 



arated variable u ■ 
tern of auxiliary classical particles. It is straightforward 
to represent the factors of the wave function ^f({uj}) in 
a product form, 
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t/>(u, t) = C(t)N(u) Y[(U- Xm(t)) , 



(2) 



where the factors C(t) and j\f(u) ensure normalization 
of the wave function. The moving roots X m (t) (m = 
1,...,M) of the wave function take time-dependence 
into account. This method for treating time depen- 
dence is due to Calogero [6 j. Substituting the prod- 
uct form ansatz ([2| into the Schrodinger equation ([!]), 
a system of first-order differential equations can be ob- 
tained: iX m = $ m ({A n }) where n, m = 1 . . . M, and the 
function ^ m ({A n }) is entirely determined by the explicit 
form of the potential A(A), ^ > m ({A n }) = a(A m )A(A m ) + 
2^1i a (^m)(A m -A n ) _1 , where a(A m ) = Ylf=i&m-Zj)' 
This system of differential equations are equations of mo- 
tion for auxiliary classical particles which move according 
to the Lagrangian 
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where x m = / dX m a 1/2 (A m ) and a(X m ) U.f =1 (X m 
zj). The corresponding Hamiltonian is therefore H 



l,m,n 

This is the classical many-body integrable model which 
governs the motion of the poles of the wave function 
and motivates a new correspondence principle for char- 
acterizing quantum dynamics: the dynamics of a quan- 
tum system can be called regular or nearly integrable 
if the tra ject ories of the auxiliary particles determined 
by Eqs. ( 3|4 ) reside on KAM tori and the breakdown 
of these tori can be seen as the criterion for irregular 
quantum dynamics. In the framework of our approach 
two different cases of perturbing the integrable system 
can be identified: one (vertical) corresponds to pertur- 
bations which leave the structure of separated variables 
untouched, while the other (horizontal) mixes separated 
variables. The first case corresponds to dynamical driv- 
ing of parameters of integrable systems while the second 
case can be achieved by adding new terms to the Hamilto- 
nian. The most generic perturbation would be a mixture 
of both. In terms of the classical system, the horizontal 
perturbation can be described by attributing a "color" 
(additional index j) to the coordinates and momenta of 
the auxiliary particles. In what follows we are going to 
test our principle on the first class of perturbations. 

Example: Dicke model - As an illustration of our for- 
malism, we consider the Dicke model which represents a 
non-trivial physical example of the class of Gaudin mod- 
els. This model, described by only one separated vari- 
able, is sufficiently generic to reveal the universal features 
of Gaudin models perturbed in the vertical direction. It 
was introduced by Dicke [14 and is widely used in the 
context of interaction of light and matter in quantum op- 
tics [15]. One particular application of the Dicke Hamil- 
tonian is the description of the Bose-Einstein condensate 
in an optical cavity [T6] . 

The model treats the interaction between a single- 
mode bosonic field (b and b\ photon annihilation and 
creation operators) and an ensemble of two-level systems 
(regarded as effective spins detuned from the bosonic 
modes by A) in the rotating wave approximation. The 
quantum Hamiltonian reads 



H D =AS z +g(tfS-+bS+), 



(5) 



^mPm + Vd^n}) with single-body, two-body and three- 



(3) where S a = J2j=i ~t ^ s a collective spin operator. The 
total number of excitations M = b^b + S z + S is a con- 
served quantity. Therefore, the relative strength of the 
detuning A/g is the only free parameter in the system 
in a given sector with well defined quantum numbers M 
and S. 
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We carry out the explicit Bethe ansatz solution of this 
model for which a(A) = A and A(A) = (A - A + 2S/X)/2 
in the Supplementary Material, where we also show that 
the underlying classical many-body system governing the 
evolution of the spectral data in this model is a so-called 
Inozemtsev-Meshcheryakov model [8 . It describes the 
motion of M particles in a nonlinear external potential, 
which interact via a long-range interaction of Calogero- 
Moser type, 

M 

Him = J2[P 2 m + V ( x ™)] (6) 

m=l 



( J ( rp rp \ 2 ( qr> I np \ 2 

n rn —\ \ n ^mj * ^m) 

n^rn 

The potential takes the form V(x m ) = -^gx^ + fix^ + 
~/x 2 m + !6S 2 x~ 2 , where /? = and 7 = \(A 2 (t) - 

A(t) — 45 + M — 1). In spite of the nonlinear potential, 
this classical model is integrable. 

Here we brake the integrability by time-dependent 
driving of the detuning. The model still pertains its 
form as in Eq. (J6|, however its parameters now depend 
on time. Concretely, we consider the following setup: 
At t = the system is prepared in its ground state at 
A = Ao. Then we evaluate numerically its time evolu- 
tion under the periodic detuning A(t) = A cos(cj£). We 
solve the dynamical equations by using a Runge-Kutta 
integration scheme. The rapidities can be obtained from 
the coefficients of the wave functions by finding the roots 
of symmetric polynomials. For these illustrative pur- 
poses, we choose a small number of excitations, M = 4, 
5 = 6 and a strong amplitude of the detuning Ao/g = 5, 
such that the bosonic modes are highly occupied initially, 
A/5 = (tfb) ~ 3.2, and the population of excited spins is 
small. The high driving amplitude causes strong dynami- 
cal redistribution of excitations between bosonic and spin 
degrees of freedom. If the driving frequency is sufficiently 
slow and non-resonant, dynamics remain almost adia- 
batic and observables are expected to exhibit periodic 
oscillations along the instantaneous ground state values. 

In Fig. 1(a) such an example is shown: At frequency 
uj = 3.57g/ft, there are regular oscillations of the bo- 
son populations N^(t) between 3.2 and 0.2. The rapidi- 
ties, which correspond to the position variables of the 
classical model (J6|, are monitored stroboscopically af- 
ter each cycle (i.e. at time t p = 27rp/co>, p = 0, . . . , P, 
where P = 4000 in the present case) by collapsing them 
onto a single complex plane. Fig. 1(c) shows that in 
this non-resonant case the rapidities cluster on circles 
located around the ground state positions. These cir- 
cles indicate the existence of stable KAM-tori in the 16- 
dimensional phase space of the classical system and ac- 
cording to our principle of correspondence between dy- 
namics of the quantum and the auxiliary classical sys- 
tem, we can classify such behavior as nearly integrable. 
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FIG. 1: Dynamics of the Dicke model driven non-resonantly 
with amplitude Ao/g = 5 and a frequency cj = 3.57g/h, 5 = 6 
and M = 4. (a) The boson occupation number Nb monitored 
over some interval of time, (b) the weights of eigenstates |7| 
Col and (c) the stroboscopic maps of all rapidities A m , m = 
1,. . . ,M after 4000 cycles. 
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FIG. 2: The Dicke model driven near- resonantly with ampli- 
tude Ao/g = 5 and frequency uj — 3.68g/h. For explanations 
see caption of Fig. 1. 
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FIG. 3: The Dicke model driven with amplitude Ao/g = 5 
and frequency uj = 3.75g/h. For explanations see caption of 
Fig. 1. 

In order to characterize the statistical properties of this 
system we measure the distribution of states averaged 
over all cycles 

c « = ^Ek^p)i«>i 2 . (7) 

p 

where \ip(t p )) is a state of the driven system at t = t p 
and \a) are the eigenstates of the Hamiltonian H(t p ) af- 
ter p cycles. In Fig. 1(b), we found their distribution 
to decay rapidly, which is expected in this nearly adia- 
batically driven system. We compare this distribution 
to a Boltzmann distribution, c a = e~P Ea /Z, with the 
same average energy. It turns out that the Boltzmann 
distribution cannot describe the weights (Fig. 1(b)). 

In Fig. 2 we consider a slight increase of the frequency 
with respect to non-resonant case to uj = 3.6Sg/H. The 
boson occupation, which starts to exhibit an additional 
beating frequency (Fig 2(a)), suggests that a resonance 
is approached in the quantum model. Interestingly, this 
comes along with a scattering of the rapidities on the 



collapsed 2-dimensional stroboscopic maps (Fig. 2(c)). 
From the point of view of the auxiliary classical system, 
these dynamics rather strongly deviate from the inte- 
grable limit. It has to be noted that despite the relatively 
dense exploration of the phase space, we could not find 
an indication of truly chaotic behavior. Nevertheless, 
and despite the small number of degrees of freedom in 
the system, this leads to a state distribution remarkably 
close to Boltzmann distribution (Fig 2(b)). 

Further increasing the driving frequency to uj = 
3.75g/h as in Fig. 3, leads to strongly beating dynamics 
of boson occupancies. This resonance of the quantum 
model leads to a new structured pattern in the strobo- 
scopic map of the classical variables. This hints that 
there are new quasi-periodic orbits emerging, which re- 
side on a topological structure different from the one of 
the near-adiabatic case. The state distribution in Fig. 
3(c) shows that there is a large amount of energy pumped 
into the system and the weights deviate considerably 
from the Boltzmann distribution. 

Conclusions and outlook - In summary, we derived 
an correspondence between a time-dependent quantum 
model and an auxiliary classical system. The strength 
of this approach is illustrated by an example of a driven 
Dicke model with a frequency tuned from a non-resonant 
to a resonant value. The emerging dynamics can be inter- 
preted in terms of the classical underlying system, whose 
trajectories lie on different KAM tori in non-resonant and 
resonant cases. At the point where one torus is deformed 
into the other, irregularity in the classical dynamics is 
most pronounced and time-averages of quantum observ- 
ables approach thermal equilibrium. 

One perspective of our method, which is not restricted 
to low-energy effective descriptions, would be an applica- 
tion to the high-energy part of quantum chromodynam- 
ics (QCD), which is described by the integrable quan- 
tum spin chain with complex spin [18]. The separa- 
tion of variables for this model has been implemented 
in Ref. [19] • Dynamical perturbation from integrability 
may therefore shed a new light on phenomena of confine- 
ment. Moreover, integrable systems found in the context 
of AdS/CFT correspondence, can provide a basis to go 
beyond the limits where the duality is established. 
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Supplementary Material 



si. SUMMARY 

A single Bethe state is not sufficient to fully describe 
the problem of the time-dependent Schrdinger equation 
of a general Bethe ansatz solvable model. Here we 
demonstrate that a sufficient basis describing dynam- 
ics of Gaudin models can be given by using quantum 
version of separation of variables, also called functional 
Bethe ansatz (FBA). We use this technique to formulate 
a theory for dynamically perturbed Bethe ansatz solvable 
models. As a result, we get a Lagrangian of a classical 
system of equations of motion which governs the motion 
of rapidities. For sake of concreteness, we focus mostly on 
Gaudin- type models, but the approach we present here 
is rather general and we show how it may be extended to 
arbitrary integrable models. 

In the first section of the Supplementary Informa- 
tion we summarize the construction the algebraic Bethe 
ansatz (ABA) and the FBA. On the basis of these tools 
we show a detailed derivation of our main results in the 
next sections. Then, we illustrate the principle of the 
method in a detailed discussion of the Dicke model. 



S2. SEPARATION OF VARIABLES IN 
QUANTUM INTEGRABLE SYSTEM 

A. Algebraic Bethe ansatz 

We start with a general construction of the ABA and 
then focus on its reduction to Gaudin models [SI . The 
ABA for the s/(2)-related models (the class of models 
where the quantum space of individual particles is a rep- 
resentation space of the si (2) algebra) starts from the 
definition of the L-operator 



The generic .R-matrix can be expanded as a regular 
series of the parameter 77, 



f A(X) B(X)\ 
L{A > ~ { (7(A) D(X) ) 



(S2) 



where the elements are operators on the quantum space 
which depend on a complex parameter A, called a 
spectral parameter or rapidity. To ensure integrability, 
the L(A)-operator is imposed to satisfy the the Yang- 
Baxter equation, 

R 12 (X - n)L 1 (\)L 2 (n) = L 2 ( f i)L 1 (X)R 12 (X - /x) , (S3) 

where the L-operators L\ and L 2 act on copies of Hi 
and % 2l and the quantum R- matrix acts on C 2 (g) C 2 . 
Both, i?-matrix and L-operators, can depend on several 
types of parameters: the spectral parameter A, inhomo- 
geneity parameters Zj (where j belongs to some countable 
set) and a parameter 77, which controls the deviation from 
the quantum group structure. This implicit dependence 
on A, {zj} and 77 is always assumed, even though only the 
spectral parameter A is kept explicitly. 



i?(A)=/ 4 x4+77r(A)+ 0(77 2 ), 



(S4) 



where the matrix r(A) is sometimes called quasiclassical 
r-matrix and 24 x 4 is the identity operator on C 2 0C 2 . 
The dependence of r on {zj} is inherited from the corre- 
sponding i?-matrix. In this limit, the Yang-Baxter equa- 
tion takes the following form: 

[L(A)0/,J0L(/i)] 

+ [r(A - /i), L(A) <8> / + / ® L(fj)] = . (S5) 

In the case of Gaudin models one assumes the following 
form for the r-matrix: 



r(A) = 



//(A) 

g(X) 

g(X) 

V /(A) 



(S6) 



Here we specialize on the case of rational functions 
/(A) = g(X) = 1/A. This leads, in particular, to spin 
models with SU(2) symmetry. Trigonometric or elliptic 
functions describe models with lower symmetry. Substi- 
tuting the matrix form of the L-operator into the Yang- 
Baxter equation one can derive an algebra of A,B,C,D- 
operators [S2 . This algebra plays a fundamental role in 
the discussion of the method of separation of variables. 
To first order in 77, and using the represent aion (S6), the 
expansion of the algebraic commutation relations (S5) 
lead to the Gaudin algebra, 

[A(X),B(»)} = f(X-(j,)B(X)-g(X-fi)B((i), 

[A{\),C(jm)] = -/(A- M )C(A) + <7(A-/x)C(/x), 

[B(X), <7( M )] = 2(/(A-m)A(A)- 5 (A-/x)A( M )), 

[B{\),Bfa)] = 0, 

[S(A),C(A)] = A'(X) - D'(X) , (S7) 

and D(X) — —A(X). The generating function of the in- 
tegrals of motion Hi up the first order in r] is defined 
as 

T(A) = ^Tr[L 2 (A)] (S8) 

= \ (A 2 (X) + D 2 (X) + B(X)C(X) + C(X)B(X)) . 

For sake of concreteness we use a realization of opera- 
tors A(A), 5(A), C(A), D(X) in terms of generators of the 
compact Lie algebra, 



N 



A W = -^/(X-z^S* 

3=0 
N 

B(X) = 5>(A-^)S7, 
3=0 

N 

C(X) = J29^-Zj)S+, 

j=0 

D(X) = -A(X) . 



(S9) 
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We explicitly restore the dependence on the inhomogene- 
ity parameters {zj} and ej. The total number of spins 
is N + 1, however, only N independent inhomogeneity 
parameters have to be considered; without restriction of 
generality one can set zq = 0. The operators are the 
generators of the Lie algebra su(2), 

[Sj, S%] = iSj5jk, [Sj, S%] = iSjdjk, [Sj, Si] = iS^5jk • 

(S13) 



We denote Sf = S] ± iS?. Substituting Eqs. p9\ into 



Eq. (S8) we obtain for the rational case 



T(A) = const + 



E 



E 



i s (i s + 1) 



(S14) 



where /j(7j + 1) = S| is the value of the Casimir operator 
of the representation on site j. The residues of T(A) 
represent Hamiltonians 



^ = E 



Sj • S/;; 

2fc 
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Particular examples are the central spin models with 
j = and 1/zj being the potential between spins, or 
Dicke models, describing interactions between bosons 
and spins, which can be obtained in the limit Iq —> oo. 

Bethe states, i.e. eigenstates of conserved quantities, 
are obtained by applying the product of C-operators to 
the vacuum state |0) (defined by B(X)\0) = 0), 

M 

|^(A 1 ,...,A M )) = l[C(\ n )\0}. 

n=l 

For the rational s/(2)-related Gaudin models the rapidi- 
ties A m , m = 1, . . . , M have to fulfill the Bethe equations 



A(A m ) = 



1 



An 



(S16) 



where the function A(A m ) depends on the set of inhomo- 
geneities {zj} and can be represented as a ratio of two 
polynomials, 



A(A„ 



E 



, A m Zj 



b + cX r 



"' Q(A m )' 

(S17) 

where dots denote possible term of higher order in A. 

B. Separation of variables 

We review the method of separation of variables to the 
classical and quantum systems. More details and further 
references can be found in the review paper of Sklyanin 
[S3] or in Ref. [S4]. 



Sklyanin [S3] suggested to introduce separated vari- 
ables for the Bethe algebra as operators uq and Uj, 
j = 1, . . . , A/", acting on the quantum space such that 



C(uj) = , and (7(A) = u /X . 



(S18) 



Arguments of A,B,C and D are no more C-numbers 
but operators. For general Bethe anstatz solvable sys- 
tems this substitution can only be done by extending the 
algebra [55], but for Gaudin- type models, the C matrix 
can be rewriten as 



N 



u rijli(A 



Uj) 



A 



A 



n.L(A 



(S19) 



and the separated variables can be given explicitly. It can 
be verified that i^o, Uj commute with each other. From 
the residue at A = zj we find that the operator Sj~ can 
be rewritten in terms of the basis of symmetric functions 
of separated variables, 



= res A ^.C(A) = ^o g^ ~. ■ (S20) 

l[k=l {Zj ~ Zk) 

A set of conjugated variables is obtained by substitut- 
ing variables Uj into diagonal elements of the L-matrix, 



Vj = A(uj) . 



(S21) 



For the algebra of Gaudin operators commutation rela- 
tions become canonical: 



[vj,u k ] = S jk . 



(S22) 



Therefore, uj and ivj can be seen as pairs of position and 
momentum operators. 

It is important to note that there is one conjugated 
pair less than spins in the system (j = 1, . . . , N). The 
operator u$ is special and does not come together with a 
conjugated operator. 

From the form of the transfer matrix T(u) = 
^Tt(L 2 (u)) and the commutation relations one can de- 
fine a realization of the Vj -operators in terms of the Uj 
variables as 



d 

dUj 



A( Uj 



(S23) 



where the function A(uj) forms a spectrum of operator 
Uj. They are given by the vacuum eigenvalues of the 
operator A(X), A(X)\0) = A(A)|0 and is therefore specific 
for every model. We note, however, that the separation 
of variables does not need a well-defined pseudovacuum 
state, like in the algebraic Bethe ansatz, and the function 
A(u) can be determined by operator techniques once we 
identified the space of states on which the operators Uj , Vj 
act |S3j . 
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1. Coordinate representation 



Using (S30) one finds 



In order to express all the operators of the model in 
terms of the new variables Uj and one considers first 
a representation of the sl(2) algebra in coordinates £j 
j = 0, . . . , N of some Euclidean space and corresponding 
differential operators d/d^: 



d 



+ lj , Sf — £j , S, 



2L 



d 



(S24) 



N II fc=l ( u j - z k) 



Of, = d. 



U 



E 



p^k 



and, from Eq. (S28) 



(S32) 



The Gaudin Hamiltonian then reads 

N 



A;=0 
N 



(S25) 



= E 



^? Z k 



N 



N N 



< = EE 



2=0 k = v J 
AT N 



(Uj - Zi)(Zi - z k ) 



(%-%) 2 , 



(S33) 



2 = fe=0 



(Uj ~ Zi)(Zi - z k ) 



1 (%-%), 

(S34) 



Since in this representation the spin operator Sj~ is 
given by the coo rdina te, differentiating the log of both 
parts of the Eq. (S19) we obtain 



duo 
u 



N * 

EdUj 
Uj - A 
i=i 3 



N 



c{\) = E 



(S26) 



In particular, this implies that 

N 



d^j duo ^ 
& ~ ^o ^ 



u k - z. 



(S27) 



Using these equations one can establish the following 
identities: 



— =y 



d 



fe=0 J 



and 







%=e4' ( s2§ ) 



k=0 



d^ = 

du k u k - Zj 



du uq 



(S29) 



These equations can be used perform a change of vari- 
ables between {^-; ^-} and {uo,uj; including 
the Jacobian. Thus, taking the pole at A = Uj a nd w rit- 
ing duj = ^Z k= i(duj / d£ k )d£ k one obtains from (S26) 



rip=i - zp) 



duj 



1 



p^j 



U ° Up=l ( U 3 ~ U p) 

P^J 



(S30) 



where the identity 

N N 



EE 

2=0 j=0 



(u — Zi)(u — Zj) 



N 



E 



N 1 

V L 



(S35) 



JV AT 



2=0 J=0 



has been used. With these formulas one finds the repre- 
sentation of the spin operator 



/ n n^=o ( u j — z^ 

u d Uo - — 

V 



S: 



3 rifcLi(^ ~ u k) 



n^=o(^ - z s ) 



1 Y\k=l (Uj ~ U k ) 



(S36) 



where we used the Eqs. (S30) and (S29) and where v 



is defined in Eq. jS23k with A(u) = V^_ n . It is 



instructive to compare it with the result of deriving it 
from the pole of A(X). Indeed, 



N 



Si 



N 

E 

fe=0 

d 

duj 



A - z k ^ 
uj - z n d£ k 



\k=0 / a=m, 



A( Uj ) 



(S37) 
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Using these formulas and Eq. (S14) together with 
Eq. ( |S8] ) we establish that the value of the generat- 
ing function T(A) taken at A = Uj (remember that 
C(uj) = 0) is 



N 



i=0 



N 

E 

i=0 



W + l) 

(Uj - Zi) 2 



J J 



d 

duj 



(S38) 



Since vj can be interpreted as a canonical variable con- 
jugated to Uj, the expression for T{uj) looks exactly as a 
Hamiltonian for the free particle propagating in the back- 
ground of some monopole (vortex) distribution in the 2D 
plane. 

The explicit expression for the commuting Hamiltoni- 
ans in terms of separated variables can be obtained ei- 
ther directly, s ubsti t uting the ex plicit change of variables 
through Eqs. p3l] ) ? ( [S32] ) into ( [S25] ) ? 



Hi = 



nN / 



Uk) 



n^=o & - z k ) 



N 

E 



nN ( 
n=0 [Ui 
n^k 



Zn) 



nN / 
71=1 {Ui 



(S39) 



Using the Lagrange interpolating formula, the gener- 



ating function for the integrals of motion in terms of 
the separated variables at arbitrary A takes the following 
form 



N 



T ^ = E n* 



Q(X)a(u k ) 



k=l 



Q'(u k )(\ - Uk)a(X) 



(S40) 



where 



N 



N 



Q(A) = JJ(A-^t), Q'{u k )= \\{u k -u p ) 

2=1 

N 



P =i 

p^k 



(S41) 



A;=0 



One can now define a generating function for the eigen- 
values of the generating function via 



N 



= E 



2 = 



A - Zi 



e(A) 
a(A) 



(S42) 



where a(A) is defined in Eq. ( |S41[ ). Therefore, the sta- 
tionary Schrodinger equation with T(A) defined in (S38) 
in separated variables takes the form of a one-dimensional 
Schrodinger equation 



d 



a(uj) ( -2A(u j )—+A 2 (u j )-A'(u j ) J V(u ,...,u N ) = e(uj)V(u , . . . , u N ) , (S43) 



\du 2 j duj , 

what proves that the wavefunction can be factorized, 

N 



9(u ,...,u N ) = Y[tl>(u j ) ( S44 ) 

3=0 

where the operator uq is now included into the separated variables. Moreover the same form of a stationary Schrodinger 



equation for the generating function T(A)\P(A) = r(A)^(A) is obtained from (S40) by taking the residues at poles at 
both sides. We note that the form of this equation is the same for every Uj, j = . . . N. 

Several remarks about the method of separation of variables are at order: (i) The separated equation for variable 
Uj depends only on the variable Uj. (ii). The separated equations have the same form for all Uj. (iii) The differential 
operator has the form of Baxter Q-operator equation, (iv) Physically, separated variables can be thought as the 
eigenmodes of the spin densities. Indeed, we can write 

8 N ga N 
j=l 3 J=0 

In other words, separated variables serve as a set of normal coordinates in the space of operators. 
Now we consider a time-dependent Schrodinger equation, 

^*({A},t) = T(A)*({A},t) (S46) 
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and make again a product ansatz, 

N 



*({A»=n^<)- ( s47 ) 

3=0 



Substituting a time-depending ansatz we obtain a separated equation 

' d 2 _ _ , x d 



idMuj) = a( Uj ) I ^ - 2A(^-)^7 + A 2 (^) - ^-'( u j) ^KO- (S48) 



Now, since the equations are the same for all ix^ we denote 
the variable Uk by u and consider the following change of 
variables 



fin) 



dy 



(S49) 



where the integration should be understood in general as 
a contour integral which avoids the singularities of a(y). 
Next we implement this transformation to the wave func- 
tion (change of variables and multiplication by the fac- 
tor). The Hamiltonian can be now put into the form of 
a one-dimensional Schrodinger equation with a Hamilto- 
nian 



H(x) = d 2 x + V{x). 



(S50) 



Instead of doing this in general we will postpone this to 
explicit examples below. We only remark that the form 
of the potential is rather special and allows for super- 
symmetrization, meaning that the function A(u) can be 
considered as a derivative o f the superpotential W(u), 
and the Hamiltonian in Eq (S48) is the bosonic part of 



the Hamiltonian of supersymmetric quantum mechanics. 

To take a time-dependence into account we will follow 
a procedure of moving zeroes as described in the book 
of Calogero [S6]. In this approach a wave function is 
assumed to have the following product form, 



A I 



(S51) 



where X a (t) is a pole of the wavefunction. For this prod- 
uct form one obtains 



d_ 
du 

d_ 

dt 

du 2 



Al 



n=l 

M 

2/j(u : t) = 2/j(u : t)J2( U - Kit))' 1 ("A n (t)) 
n=l 

M 

tl>{u,t) = 2^,t)^(w-A n (t))- 1 



A I 



]T (\ n (t) - x m (t)y 



(S52) 



Substituting this into Eq. ( S48 ) and equating the residues 
at the poles we obtains 



iX n (t) = 2a(X n ) 



1 



rn^n 



An — X r , 



2a(A n )A(A n ). (S53) 



At equilibrium, A n = 0, we recover the Bethe equation 
in the form of Eq. (S16). 



Now we are goinge to derive that equations ( S53 ) can 



be considered as equations of motion for certain La- 
grangean dynamical s ystem . Indeed if we denote the 
right hand side of Eq. ( S53| ) by $ n ({A}) and differentiate 
the equation —iX n (t) = <£ n ({A}) we obtain 



M •VK 

— * / J ~W\ *n 

^ dX m 

m=l 



A I 



(S54) 



On the other hand if we consider the Lagrangean for M 
particles 



A I 



A I 



£=£(*n) 2 -£(*n({A})) 2 



(S55) 



n=l 



n=l 



regarding the set {X n (t)} as generalized coordinates, on 
can derive the equations of motion, 



A n 



1 dC 

2dX n 



M 



E*. 



OK 



(S56) 



We note that equations, (S54) and (S56) hold only if 



d\ r , 



0. 



(S57) 



d is an exterior derivative. For the gaudin models, one 
can show that <I> n = d n (j)({Xk}) with some function 
0({A&}), therefore the 1-form is always closed. If one 
makes the change of variables ( |S49| ) from variable A n to 
variables x n the equation of motions acquire the form 
— ix n = d Xn (j) Such, one establishes a M-particle La- 
grangean system whose Hamiltonian is given by 



A I 



H = Y^vl + V M {{x n }). 



(S58) 



n=l 
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where p n = x n . Substituting directly the function 
<E>({A}) into the Lagrangean gives for the potential, 
in terms of {A n } (explicitly obtained by expanding 
V„i -1-7, ./„)). 



A I 



V M ({\}) = 4^a„A2+8 Yl 



CL n A n 



n—l \m^n J 

where a n = a(A n ) and A n = A(A n ). Note that after 
expressing this in terms of A n 's we have to make a a 
change of variables to x — J a~ 1 / 2 (\)d\ and express Vm 
in terms of them. For generic functions a n (A), A n (A) we 
may get rather complicated potential. We will show an 
explicit form of the potential Vm for the Dicke model 
below which gives an integrable classical model. 



S3. SOLUTION OF THE DICKE MODEL AT 
EQUILIBRIUM AND WITH TIME-DEPENDENT 
DETUNING 

A. Equilibrium BA 

In this section we are going to discuss the details of the 
mapping dynamics of the quantum Dicke Hamiltonian 



H D = AS Z +g(tfS~ +bS+), 



(S60) 



to the classical equations of motion. The total spin is 
Y^a=i(S a ) 2 = S(S + 1) and the conserved quantity is 
M = tfb + S z + S. It can be shown [S7], that this model 
belongs to the Gaudin algebra with 



C(A) = 6 f 



St 



A 



(S61) 

|0) = |0)6|S,-S). (S62) 
At equilibrium the Bethe equations read 

M 



25 

A 77, 



A„ + A- 



m=l 



(S63) 



and define the eigenenergies 

M 

X 



M M 

E SM = -S(A + 2 J2 t- ) = A(M -S)-J2 V 



=i 



(S64) 



B. Time-dependent perturbation 



A natural way to drive the Dicke model out of equi- 
librium is to apply some time-dependent perturbation. 



For time-dependent detuning we obtain a solution of 
the Schrodinger equation from a Bethe wavefunction in 
which static spectral parameters are replaced by time- 
dependent ones and then derive the underlying classical 
system of equations of motion. We make the ansatz 

M 

|tt(t)> = exp[ie(t)] JJ C(A n (t))|0> , (S65) 

n=l 

where 

rt 



e(t) 

6(/i,r) - ^ 2 (/i) 



/ res /i ^ S(M,^)(ir , 
Jo 



dx(/i) 



M 



X(fl) - x(\ n (T)) 

M - A„(r) 



e(t) 



,t M 

S / (A + 2V[A„(r)]- 1 )dr. 

J n=l 



(S66) 



e(t) is the time-dependent energy. The function x{fi) is 
defined by the vacuum state. The spectral parameters 
are now subject to the following set of equations 



,. A n (t) 
! A„(i) 



A/ 



1 



; A n (i) -x m (t) 



(S67) 



In what follows we show that the corresponding classi- 
cal auxiliary system is the Inosemtsev model. After the 
change of variables ( S49 ) 

\ n = x 2 l , x n ^x n /2, (S68) 

the equation for the dynamical BA equation takes the 
following form: 

r n = ax\ + bx n + — + 2x n ^ ^ 2 ^ ^ 2 . (S69) 



The constants are a = 1/8, 6 = —A/2, c = —46'. Con- 
sidering now the equation 

^2 



(2ac + 6 2 + 2a(M - 1)) x n - 4abx^ - 3a 2 xl 



M 

■ E 



one recovers the potential of the classical model: 

M 2 i 



(S70) 



(S71) 



n=l 
M 



^(Kl) = E ( 



1 



1 



l M 



16S 2 



n=l 
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The parameters are given by 

P = -^jT> j=\(A 2 (t)-A(t)-4S + M-l). 

(S72) 



This model is known as BC-type Inozemtsev model [S8 . 
It is integrable on the classical level for time-independent 
parameters. We note that this is a complexified version of 
the classically integrable Inozemtsev (or BC-type) model. 
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